Intro 1st paragraph:

Large scale envrionmental changes may drive changes towards different vegetation types. Biomes with alternative stable states that could be tipped into different vegetation types with small changes in environment create uncertainty in the future vegetation response, and fuel fears of forest instability (Staver et al. 2011). Globally, savannas and forests are often cited to be alternative stable states. However, there is widespread debate over whether alternative savanna and forest states exist or are an artifact of satellite processing (CITE Modis paper, STAVER response), and high uncertainty in predicting their sensitivity to environmental changes (Higgins and Scheiter 2012). Therefore, deterimining if savanna and forests existed as alternative stable states in the recent past and evaulating their past responses to environmental changes can provide key insights into the impact of future changes on ecosystems with alternative states. To explore the potential responses of alternative stable states to environmental changes, we use historical (pre-European Agricultural Settlement) and modern vegetation data from the North American Midwest to ask whether there is evidence for savanna and forest alternative stable states in the temperate zone and determine if these conditions have been resilient to large scale fire suppression, changes in climate. 

In the North American temperate zone, much of the savanna-forest boundary has transitioned into a closed forest state (non-agricultural landscape) (Nowaki and Abrams 208). This shift has occured over a time period of large scale fire supppression at least since the 1930's (Nowacki & Abrams 2008). We estimate pre-European agricultural settlement stem density across the Upper Midwest using historical survey data collected in the 19th century (Public Land Survey (PLS)). We evaluate whether there is evidence of alternative stable states in the past by determining if the past landscape had a bimodal distribution of tree density that cannot be accounted for using climate or soils. We then assess how persistant the hypothesized alternative stable states are on the landscape by comparing the historic relationships between density, soils, and climate to the modern relationship of tree denisty to soils and climate. Finally, we apply the past and modern relationships with soils and climate to assess the resilience to projected future climatic changes. We hypothesized that  PLS distribution of tree density was bimodal and not well explained by environmental factors due to latent intact vegetation-disturbance feedbacks along the savanna-forest boundary before settlement, but fragmentation and fire suppression has stabilized modern vegetation and produced a modern boundary that is not bimodal.

Understanding forest response to environmental changes is a key challenge to predicting and planning for the future. Vegetation response may be a deteriministic, smooth response to climate in some places, but non-linear and non-deterministic vegetation-climate relationships can result in divergent forest responses. Globally, transistions from savanna to forests are sharp, and often not deterimined soley by climate. Rather, savannas and forests may be “alternative stable states” that coexist within the same climate due to vegetation structure-disturbance feedbacks (Staver et al.2011, CITE others). Savanna and forest alternative states are largely observed in the tropics, fueling fears that the biomes are potentially unstable and that environmental changes and increased disturbances could result in a a collapse of closed forests to the low tree density savanna state (Staver et al. CITE others). The potential collapse of forests that can switch to their alternative savanna stable state would have enormous consequences for terrestrial carbon storage and biodiversity. Similar vegetation structures (savannas and forests) are also observed in the temperate zone, for example, historically, in the North American midwest, the Humid Pampas region in South America, and (temperate steppe in Asia), suggesting a potential that vegetation structure-disturbance feedbacks could have lead to alternative stable states in these regions as well. However, large scale environmental changes (land use, fire suppression, climate changes) have already occurred in the temperate zones, making these regions ideal to study the resilience potential savanna/forest instability to past environmental changes. [paragraph purpose: savannas and forests as alt. stable states, temperate tropics potentially as well, ??uncertainty in envtl. response, we could study this in the temperate zone]

While the climates of temperate and tropical forest/savanna systems differ, the proposed mechanisms for savanna and forest vegetation in both regions converge on disturbance regimes (however, this is well studied in the tropics, where disturbance regimes are intact). While climate does not deterministically predict the vegetation state overall in these regions, precipitation amount, seasonality, and soil factors can all play a role (Lehmann et al 2011, Staver et al. 2011, Danz et al. 2011). Particularly at a large regional scale, deterministic high and low precipitation creates forests and open savs./grasslands. respectively (Danz et al. 2011, Staver et al. 2011). Additionally, hetergeneity in soil factors that affect hydrology and overall climate may be important at smaller spatial scales (Danz 2011). Soil and topography are hypothesized to act as fire breaks, creating upland forested regions and lowland, flat prairies in the temperate zone (grimm, nowacki & abrams). Vegetation structure-disturbance feedbacks are important intermediate precipitation levels in seasonally dry tropical regions (Staver et al. 2011). In these intemediate climates, fire is associated with open savanna ecosystems, where high grassy fuel loads in open systems promote fire spread, which prevents understory trees from establishing. However at these same climates, closed forests may be present but light limition in the understory prevents fuel buildup, suppressing the likelihood of fire (@hoffmann_2012). These vegetation-disturbance feedbacks drive will drive regions of intermediate tree cover either towards a closed forest state in the absence of disturbance, or towards an open savanna state in the presence of disturbance. Thus, savanna and forest stable states in these regions often results in a hallmark bimodal distribution in tree cover that is not deterministically predicted by climate or environmental factors. The instability and stochastic nature of alternative stable states and disturbance feedbacks make it difficult to predict the outcome of these vegetation processes into the future. [Paragraph Purpose: Alternative stable state mechanism and background, temperate and tropical hyp. mechanisms are extremely similar]

In this study, we use historical vegetation data from the North American Midwest to ask whether there is evidence for savanna and forest alternative stable states in the temperate zone and determine if these conditions have been resilient to large scale fire suppression, changes in climate.

Methods:

*Public Land Survey Data:* The Public Land Survey was commissioned by the U.S. General Land Office to assess the quality of land/timber and assign land titles for Euro-american settlers. The surveyors placed corner posts at 1 mile section corners in a grid within each township. At these section corners, and at 1/4 section corners, the surveyors would record the distance and azimuth from the corner to the nearest two trees, the species of the nearest two trees, and the diameter at breast height (DBH) of those trees. After these historic records were digitized (cite Mladenoff et al. ), we estimate the denisty and basal area at each corner point using the unbiased morista density estimator with correction factors for surveyor bias (@goring_2016, CITE Cogbill). Tree density is averaged across 8km x 8km grid cells in the upper Midwest.

*Forest Inventory Analysis Data:*

Add FIA data methods from Sean/Simon/Jack. FIA estimated tree denisty is also aggregated to the same 8km raster grid resolution.

Climate Data

Mean annual temperature and Mean annual Precipitation data for the modern and historic vegetation-climate relationships are from PRISM datasets. Modern mean annual temperature and precipitation fore each 8km grid cell is calculated from the 800m 30-year Normals dataset. Historical climate is estimated as the average annual Temperature and average annual Precipiation over the 1895-1910 period. While this historical climate period does not overlap with the entirety of the European agricultural settlement era, this climatic period is the oldest for which we have reliable estimates of temperature and precipitation in the region. Precipitation seasonality was calculated from Monthly Prism datasets as follows:

SI = sum 1-12(abs(mi - P/12))/P * 100

Temperature is calculated using the coefficient of variation method TSI = sd(m1….m12)/Tavgannual *100

Soil data

Soils data are derived from statewide gSSURGO 10m raster data product (https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/home/?cid=NRCS142P2_053628). The weighted averages from the top 30cm of the soil were calculated for % sand, available water content (awc), and saturated hydraulic conductivity (ksat) soil parameters using the gSSURGO On-Demand ArcGIS toolbox (https://github.com/ncss-tech/ssurgoOnDemand.git). Weighted averages were calculated on the statewide soil raster datasets for Minnesota, Michigan, Wisconsin, Indiana, and Illinois, then weighted average rasters were mosaicked together to produce one single raster dataset. Maps of soil parameters were then aggregated to a 1km grid scale and an 8km Great Lakes St. Lawrence projected grid scale (PalEON dataset scale).

Pricipal Component Analysis

Since many climate and soils variables are often correlated, a pricipal components analysis was conducted to reduce dimensions and co-linearity in the environmental data. PCA was coducted on the scaled variables of historic MAP, historic Mean temperature, historic Precipitation Seasonality Index, Historic Tempearture Seasonality, AWC, and %sand. The modern principal component scores for climatic and soil variables were predicted using the historic principal component model. This method maintains the PCA of modern and historic environmental data on the same scale. The 1st Principal Component (PC1) was then used as a covariate when assessing bimodality.

Bimodality analysis

Criteria for bimodality uses a both a bimodality coefficient (BC) and Hartigans Diptest for unimodality. Bimodality Coefficient is calculated from size, skewness, and kurtosis of the distributions (Pfister et al. 2013). A BC greater than 0.556 suggests bimodality. Bimodality coefficient (BC) is calculated from the distribuiton of data as follows using the R package ‘modes’:

BC = \(\frac{(skew^2 + 1)}{(kurtosis + 3)*\frac{(n-1)^2}{(n-2)*(n-3)}}\)

Desnity distributions of the PLS and FIA data were estimated. We then calculated the BC and performed Hartigan’s diptest for unimodaltity for the on these smoothed density distribuitons for the PLS and FIA overall, and the BC within different bins of precipitation and the 1st Principal Component. Hartigan’s diptest has a null hypothesis of a unimodal distribution of the data (CITE). This provides one BC value for a given range of precipitation/PC and whether the distribution is significantly different from a unimodal distribution. These ranges are used to determine the climate space where savannas and closed forests coexist.

Generalized Additive Models

In addition to evaluating the BC for ranges of data, we developed several generalized additive models (GAMs) to determine if tree density can be deterministically predicted by climate and/or soils. GAMs were developed with both strictly additive terms and with flexible smooth terms. The full PLS data and FIA data were separated into two random halves that make up the testing and training datasets. GAM’s were developed with the training dataset, and prediction error was estimated using the test dataset. Model parsimony was evaluated basing on AIC values, and model fit was assessed using RMSE.

Future Forest Stability under CMIP5 projected climate scenarios

To investigate the consequences of the modern climate-vegetation relationship, and the past climate-vegetation relatiohsip would have for the stability of future midwestern forests, we utilized climate projections from CMIP5 representative concentration pathways (rcps) 2.6, 4.5, and 8.5. We obtained average monthly precipitation, total annual precipitation, average monthly temperature, and average annual temperature projections for 2070-2099 under the three RCP scenarios from the CCESM4 GCM climate downscaled climate projections (from worldclim, cite). We calculated precipitation and temperature seasonality as described above. Principal componenet scores (based on the past climate PCA) were assigned for each grid cell using the projected climate data from each RCP.

We then identified the climate space (within PC1) that was bimodal in the PLS time, and mapped where this climate space is projected to be in 2070-2100. We also identified the climate space where FIA forests were bimodal on the modern landscape, and where this climate space would be in 2070-2100.

Results:

Distribution of modern and past data

The distribution of tree density in the past is significantly bimodal across the entire upper midwest, with a prominant low tree density mode at 30 trees/ha and a high tree density mode at 205 trees/ha (Diptest p < 0.05, BC = 0.74306). While Hartigan’s diptest shows shows a significant multimodality to the Modern FIA tree denstity, the two modes both occur at intermediate and high tree denisty (147 trees/ha, and 550 trees/ha). Therefore, we do not see a signifcant bimodality with low and high tree denisty models on the modern landscape. This shift in tree density is well explained by the original vegetation. Specifically, grid cells with historically low tree density have generally increased in tree denisity, while those that were dense closed forests in the past have decreased in tree density (Linear reg., R^2 = XXXX, p < 0.05, Figure 2). In the past, grid cells with low tree denisty were comprised of mostly Oak and Hickory species. However, grid cells on the modern landscape with low species density generally have mixed species composition of mixed or single species hardwoods. Most low density FIA grid cells are not dominated by Oak or Hickory species, representing a shift in species composition that accompanied the shift in tree denisty.

PCA (supplemental material??)

Principal component 1 (PC1) explains 54% of variance and component 2 (PC2) explains 29% of the variance. Past precipitation and Past Temperature both have strong positive loadings on PC1, while temperature seasonality and precipitation seasonality index both have negative loadings fro PC1. Sand has the strongest negative loading and AWC has strongest positive loading on PC2. Plotting density/vegetatation type in Principal component space does not separate vegeatation into savanna, prairie, and forest.

Generalized Additive Models/Relationship between historic vegetation and climate/soil variables

Gam model selection for the model with the lowest AIC value yeilded the model: tree density ~ s(MAP) + s(mean annual temperatue) + s(deltaT) + s(deltaP). While AIC suggests this is the best model (Table 1), it only explains 61% of the Deviance and has high Mean squared prediction error for predicting the test data set (MSPE = 5976.366). Prediction of tree density is overestimated in low tree denisty grid cells and underestimated in high tree denisty grid cells.

The Gam model for the modern landscape explained only 29 % of the deviance and had a high mean squared prediction error for predicting the test data set (MSPE = XXXX). Prediction of tree density is again overestimated in low tree denisty grid cells and underestimated in high tree density grid cells.

Notes: perhaps we should find the places that there is bimodality and run the gams in these places to see if soils or something else might predict tree density within these bimodality regions.

Historic climate relationship

PLS tree density was significantly bimodal between xx - xx scores of PC1, which corresponded to intermediate values of temperature and precipitation (XXXX mm/year and XX - XX degC). Over 34% (212,032 km^2) of the the past landscape was bimodal (figure 4) during the PLS time.

Modern climate relationship Modern tree denisty was not significantly bimodal overall, and bimodality was detected at less than 1% of the landscape. However, if the historic climate vegetation relationship had remained intact, and only climate had shifted, we would have expected 34.5% of the modern landscape to be bistable today (Figure 4).

Discussion:

Tree density in the past was significantly bimodal, and was not accounted for by climate or soil factors, supporting our hypothesis that savannas and forests were alternative stable states in the region. Projecting the past climate-vegeation relationship onto the modern climate would drive shifts in the climate space where bistable tree distributions are possible, but continue to predict large extent of savanna and forest alternative stable states on the modern landscape. However, modern tree density is no longer bimodal on the modern landscape, suggesting an unexpected shift towards stable closed forests in the region. Neither past tree density nor modern tree density has a fully deterministic relationship with climate, but bimodality in past tree density is predominant at intermediate temperature and precipitation space. These largescale shifts from alternative savanna and forest stable states towards stable closed forests indicate that land use changes (fire suppression, logging, agricultural conversion) have had a far greater impact on midwestern woody ecosystems than changes in climate over the last ~150 years.

Land use changes that have led to fire suppression and land conversion in the region have long been hypothesized to drive “mesophication” (Nowacki and Abrams 2008), and increases in woody cover in the region. Additionally, small scale studies indicate that historic geography that would act as fire breaks (rivers and topography) are important explanatory variables for PLS vegetation (Grimm 1984). Modern land management strategies for restoring open savannas and grasslands in the region require extensive cutting and prescribed burns to limit understory tree growth (CITE). (CITE) suggest that a shift in overall fire return intervals have shifted southern midwest savanna and grasslands towards mesophytic forests and require larger and longer efforts of prescribed burns to return these ecosystems to their open states. While there has been small scale evidence that fire and disturbances were historically important for the maintence of temperate savanna-forest boundary in the US, our study is the first to document evidence for alternative savanna and forest stable states across a large region of the midwest. Additionally, the modern landscape after long term land use changes and post-settlement fire suppression has become mostly stable forests.

Evidence for alternative savanna and forest stable states in the temperate Midwest provides evidence that bistable systems are not constricted to the tropics. In the tropics, intermediate precipitation and high precipiation seasonality is the climatic domain for alternative stable states. However, our finding of temperate alternative stable states at much lower precipitaiton levels and temperatures, indicate that the constraints on growth identified in the tropics do not impose hard limits on the extent of these systems globally. Rather, we propose that alternative stable states are constrained to regions of intemediate moisture, where disturbance feedbacks occur.

Overall the shifts between the modern landscape and the pre-European settlement is marked by large scale forest stabilization, a shift from bimodal distribution in tree cover characterized by closed forests and open savannas, to mostly closed forested ecosystems. Throughtout the holocene, shifts in tree composition and tree cover have been linked to climatic variability across the region (CITE). However, this recent shift is not predicted by the observed changes in temperature, precipitation, or T/P seasonality. Rather, we propose that increased fire suppression and removal of disturbances from the landscape as likely mechanisms for overall shift towards closed forests. However, lack of direct evidence for the mechanism of these shifts makes it difficult to direct attribute the cause. Additonally, increases in atmospheric CO2 have been hypothesized to have a CO2 fertilization effect on tree growth, and could increase forest cover by increasing growth and tree recruitment (Wycoff and Bowers, CITE recruitment paper). To determine how much a CO2 fertilization effect may have contributed to increased forest cover in the Midwest, mechanistic evidence for CO2 fertilization in the region is needed.

The distribution of tree density across the region and the stabiliity of ecosystems to envrionmental changes can have large consequences for future biodiversity, carbon storage, and management decisions. Although evidence for alternative savanna and forest stable states in the midwest and throughout the globe suggest the potential for large scale forest/biome instability with changes in climate, ongoing land use changes and fires suppression have historically stabilized forests in the midwest. While increased stability of terrestrial ecosytems is a positive effects of the “mesophication” and overall conversion of savannas and grasslands to forests, it may come at the large cost of species biodiversity loss, increased vegetation homogenity, and loss of habitat. Assuming current land cover strategies continue, the modern forests are not likely to revert back to savanna in response of climate. However, increased likelihood of fire disturbances and further land use changes that may accompany changes in climate could alter the stability of the region. Being able to identify forest stability, and vegetation responses to different types of environmental changes will be important to understanding the trajectory of future forest and savanna biomes.

Discussion outline:

Figures:

Figure 1:Public land survey tree density is bimodal and FIA tree density is not:

Figure 2: Change in tree density between the modern and historic vegetation depends on original vegetation density.

Figure 3:

Using the climate space with reduced dimensionality (PC1) and the vegetation classificaitons from rheumtella, we can map out places where tree density was bimodal in the past:

Figure 4a:

Figure 4b:

Supplemental Figures:

Additionally, the biplot for the Principal Components Analysis:

Climate space where the bimodality occurs: (3d plots):

Table 1: PLS density GAM models for stable sites 

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## PLSdensity ~ s(pasttmean) + s(MAP1910)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  191.284      2.843   67.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                edf Ref.df     F p-value    
## s(pasttmean) 8.433  8.903 17.55  <2e-16 ***
## s(MAP1910)   7.169  8.237 33.81  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.344   Deviance explained = 35.2%
## GCV =  10094  Scale est. = 9958      n = 1232
  model formula df AIC MSPE dev.expl
PLS.gam15 PLS.gam15 PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 42.94 14186 5919 63.93
PLS.gam14 PLS.gam14 s(sandpct) + s(awc) 39.43 14203 5898 63.21
PLS.gam13 PLS.gam13 PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 26.47 14365 6618 57.16
PLS.gam12 PLS.gam12 s(sandpct) + s(awc) 24.05 14435 7124 54.48
PLS.gam14L PLS.gam14L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 6 14616 8111 45.7
PLS.gam15L PLS.gam15L s(sandpct) + s(awc) 6 14616 8111 45.7
PLS.gam12L PLS.gam12L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 14622 8165 45.36
PLS.gam11 PLS.gam11 s(sandpct) + s(awc) 27.93 14803 9093 39.03
PLS.gam10 PLS.gam10 PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 20.08 14817 9304 37.54
PLS.gam9 PLS.gam9 s(sandpct) + s(awc) 25.66 14853 9542 36.25
PLS.gam7 PLS.gam7 PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 17.6 14857 9668 35.22
PLS.gam1 PLS.gam1 s(sandpct) + s(awc) 9.717 14997 10605 26.44
PLS.gam13L PLS.gam13L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 15005 10986 25.41
PLS.gam10L PLS.gam10L s(sandpct) + s(awc) 5 15017 10995 24.68
PLS.gam11L PLS.gam11L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 6 15018 10967 24.77
PLS.gam8 PLS.gam8 s(sandpct) + s(awc) 9.296 15054 11186 22.94
PLS.gam9L PLS.gam9L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 15080 11685 20.74
PLS.gam7L PLS.gam7L s(sandpct) + s(awc) 4 15098 11878 19.42
PLS.gam2 PLS.gam2 PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 10.68 15101 12069 20.13
PLS.gam5 PLS.gam5 s(sandpct) + s(awc) 5.492 15126 12040 17.76
PLS.gam3 PLS.gam3 PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 9.658 15127 12161 18.25
PLS.gam4 PLS.gam4 s(sandpct) + s(awc) 10.15 15158 12365 16.29
PLS.gam8L PLS.gam8L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 4 15175 12246 14.24
PLS.gam5L PLS.gam5L s(sandpct) + s(awc) 3 15180 12384 13.76
PLS.gam2L PLS.gam2L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 15180 12832 13.73
PLS.gam6 PLS.gam6 s(sandpct) + s(awc) 10.5 15237 13271 10.78
PLS.gam4L PLS.gam4L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 15308 14076 4.33
PLS.gam6L PLS.gam6L s(sandpct) + s(awc) 3 15314 14027 3.832
PLS.gam1L PLS.gam1L PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 15318 14176 3.572
PLS.gam3L PLS.gam3L s(sandpct) + s(awc) 3 15335 14007 2.196

Table 1A: PLS denisty GAM models with log transformed PLSdensity as as the dependant variable 

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(PLSdensity) ~ s(pasttmean) + s(MAP1910)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.88276    0.02488   196.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                edf Ref.df     F p-value    
## s(pasttmean) 4.055  5.058 34.79  <2e-16 ***
## s(MAP1910)   8.365  8.878 58.83  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.362   Deviance explained = 36.9%
## GCV = 0.77084  Scale est. = 0.76244   n = 1232
Table continues below
  model formula df AIC MSPE
PLS.log15 PLS.log15 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 43.08 2656 0.5262
PLS.log14 PLS.log14 s(sandpct) + s(awc) 38.88 2677 0.5305
PLS.log13 PLS.log13 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 27.78 2725 0.5611
PLS.log12 PLS.log12 s(sandpct) + s(awc) 22.97 2884 0.6378
PLS.log11 PLS.log11 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 29.54 2996 0.7545
PLS.log10 PLS.log10 s(sandpct) + s(awc) 24.9 3011 0.7629
PLS.log12L PLS.log12L log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 3065 0.7532
PLS.log14L PLS.log14L s(sandpct) + s(awc) 6 3066 0.7541
PLS.log15L PLS.log15L log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 6 3066 0.7541
PLS.log9 PLS.log9 s(sandpct) + s(awc) 23.71 3095 0.7891
PLS.log7 PLS.log7 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 14.42 3177 0.8383
PLS.log10L PLS.log10L s(sandpct) + s(awc) 5 3274 0.9129
PLS.log11L PLS.log11L log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 6 3276 0.9128
PLS.log1 PLS.log1 s(sandpct) + s(awc) 10.12 3341 0.9737
PLS.log9L PLS.log9L log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 3374 0.9777
PLS.log13L PLS.log13L s(sandpct) + s(awc) 5 3374 0.9901
PLS.log8 PLS.log8 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 16.11 3400 1.045
PLS.log7L PLS.log7L s(sandpct) + s(awc) 4 3434 1.023
PLS.log5 PLS.log5 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 10.6 3447 1.09
PLS.log8L PLS.log8L s(sandpct) + s(awc) 4 3502 1.133
PLS.log3 PLS.log3 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 9.806 3506 1.095
PLS.log5L PLS.log5L s(sandpct) + s(awc) 3 3509 1.139
PLS.log2 PLS.log2 log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 9.994 3589 1.184
PLS.log4 PLS.log4 s(sandpct) + s(awc) 10.55 3603 1.215
PLS.log2L PLS.log2L log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 3615 1.228
PLS.log6 PLS.log6 s(sandpct) + s(awc) 9.716 3638 1.229
PLS.log6L PLS.log6L log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 3667 1.274
PLS.log3L PLS.log3L s(sandpct) + s(awc) 3 3677 1.278
PLS.log4L PLS.log4L log(PLSdensity) ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 3687 1.289
PLS.log1L PLS.log1L s(sandpct) + s(awc) 3 3717 1.33
  dev.expl
PLS.log15 60.54
PLS.log14 59.61
PLS.log13 57.21
PLS.log12 50.94
PLS.log11 46.84
PLS.log10 45.81
PLS.log12L 41.5
PLS.log14L 41.54
PLS.log15L 41.54
PLS.log9 41.86
PLS.log7 36.89
PLS.log10L 30.66
PLS.log11L 30.66
PLS.log1 27.44
PLS.log9L 24.86
PLS.log13L 24.85
PLS.log8 24.6
PLS.log7L 20.98
PLS.log5 20.94
PLS.log8L 16.48
PLS.log3 16.98
PLS.log5L 15.84
PLS.log2 11.24
PLS.log4 10.27
PLS.log2L 8.292
PLS.log6 7.589
PLS.log6L 4.355
PLS.log3L 3.558
PLS.log4L 2.798
PLS.log1L 0.3809

Table 2: FIA density GAM models

  model formula df AIC MSPE
FIA.gam15 FIA.gam15 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 31.27 19006 27066
FIA.gam11 FIA.gam11 s(sandpct) + s(awc) 21.87 19298 14484
FIA.gam10 FIA.gam10 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 14.32 19321 14229
FIA.gam14 FIA.gam14 s(sandpct) + s(awc) 27.44 19331 26678
FIA.gam10L FIA.gam10L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 5 19346 16207
FIA.gam12 FIA.gam12 s(sandpct) + s(awc) 15.18 19347 25459
FIA.gam11L FIA.gam11L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 6 19348 16240
FIA.gam16 FIA.gam16 s(sandpct) + s(awc) 11.96 19414 14248
FIA.gam14L FIA.gam14L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 6 19421 23095
FIA.gam15L FIA.gam15L s(sandpct) + s(awc) 6 19421 23095
FIA.gam12L FIA.gam12L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 5 19437 25197
FIA.gam16L FIA.gam16L s(sandpct) + s(awc) 4 19444 16088
FIA.gam3 FIA.gam3 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 10.07 19457 31573
FIA.gam13 FIA.gam13 s(sandpct) + s(awc) 22.38 19462 20892
FIA.gam8 FIA.gam8 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 12.03 19490 11629
FIA.gam5 FIA.gam5 s(sandpct) + s(awc) 8.802 19534 11648
FIA.gam8L FIA.gam8L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 4 19542 10886
FIA.gam5L FIA.gam5L s(sandpct) + s(awc) 3 19547 10772
FIA.gam3L FIA.gam3L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 3 19566 29092
FIA.gam9 FIA.gam9 s(sandpct) + s(awc) 20.37 19645 14473
FIA.gam13L FIA.gam13L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 5 19653 14621
FIA.gam9L FIA.gam9L s(sandpct) + s(awc) 5 19691 15974
FIA.gam7L FIA.gam7L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 4 19703 16435
FIA.gam1 FIA.gam1 s(sandpct) + s(awc) 8.143 19765 13264
FIA.gam1L FIA.gam1L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 3 19805 18718
FIA.gam6 FIA.gam6 s(sandpct) + s(awc) 10.29 19873 12076
FIA.gam4 FIA.gam4 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 7.758 19888 12145
FIA.gam4L FIA.gam4L s(sandpct) + s(awc) 3 19889 12690
FIA.gam2 FIA.gam2 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 9.375 19897 11659
FIA.gam6L FIA.gam6L s(sandpct) + s(awc) 3 19912 11451
FIA.gam2L FIA.gam2L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 3 19913 11741
FIA.gam7 FIA.gam7 s(sandpct) + s(awc) 16.35 NA 14956

Table 2A: FIA density GAM models modeled with log transformation of FIA density

  model formula AIC MSPE dev.expl
FIA.gam15 FIA.gam15 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19006 Inf 33.74
FIA.gam11 FIA.gam11 s(sandpct) + s(awc) 19298 Inf 19.15
FIA.gam10 FIA.gam10 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19321 5.034e+288 17.15
FIA.gam14 FIA.gam14 s(sandpct) + s(awc) 19331 Inf 32.78
FIA.gam10L FIA.gam10L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19346 Inf 14.83
FIA.gam12 FIA.gam12 s(sandpct) + s(awc) 19347 Inf 31.06
FIA.gam11L FIA.gam11L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19348 Inf 14.84
FIA.gam16 FIA.gam16 s(sandpct) + s(awc) 19414 Inf 11.83
FIA.gam14L FIA.gam14L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19421 Inf 26.94
FIA.gam15L FIA.gam15L s(sandpct) + s(awc) 19421 Inf 26.94
FIA.gam12L FIA.gam12L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19437 Inf 26.1
FIA.gam16L FIA.gam16L s(sandpct) + s(awc) 19444 8.606e+291 9.211
FIA.gam3 FIA.gam3 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19457 Inf 25.65
FIA.gam13 FIA.gam13 s(sandpct) + s(awc) 19462 Inf 26.56
FIA.gam8 FIA.gam8 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19490 Inf 7.458
FIA.gam5 FIA.gam5 s(sandpct) + s(awc) 19534 Inf 4.385
FIA.gam8L FIA.gam8L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19542 2.339e+273 3.293
FIA.gam5L FIA.gam5L s(sandpct) + s(awc) 19547 4.656e+233 2.859
FIA.gam3L FIA.gam3L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19566 Inf 19.62
FIA.gam9 FIA.gam9 s(sandpct) + s(awc) 19645 Inf 17.39
FIA.gam13L FIA.gam13L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19653 Inf 15.36
FIA.gam9L FIA.gam9L s(sandpct) + s(awc) 19691 Inf 13.27
FIA.gam7L FIA.gam7L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19703 Inf 12.49
FIA.gam1 FIA.gam1 s(sandpct) + s(awc) 19765 1.12e+207 9.527
FIA.gam1L FIA.gam1L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19805 Inf 6.572
FIA.gam6 FIA.gam6 s(sandpct) + s(awc) 19873 2.046e+286 3.412
FIA.gam4 FIA.gam4 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19888 1.059e+178 2.159
FIA.gam4L FIA.gam4L s(sandpct) + s(awc) 19889 6.106e+189 1.543
FIA.gam2 FIA.gam2 FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19897 1.182e+210 1.819
FIA.gam6L FIA.gam6L s(sandpct) + s(awc) 19912 4.357e+166 0.06073
FIA.gam2L FIA.gam2L FIAdensity ~ s(MAP2011) + s(modtmean) + s(moddeltaT) + s(moderndeltaP) + 19913 8.102e+165 0.009375
FIA.gam7 FIA.gam7 s(sandpct) + s(awc) NA 7.22e+296 17.1

Predicting Tree density (GAMs)

The model with the lowest AIC value is represented by the formula: PLSdensity ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + s(sandpct) + s(awc)

Supplemental PLS density model fit figures:

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -72.75  107.10  199.40  187.10  267.50  437.90
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.664  78.470 174.600 175.500 261.200 757.600
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   85.01  185.50  284.50  278.60  368.70  495.90     114
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   85.01  185.50  284.50  278.60  368.70  495.90     114

The predicted vs. observed plot shows overestimation of PLS tree density in Low density areas and underestimation of PLS tree density in High density areas. Below we map out the predictions and prediction errors in space. Model fit has improved for the stable portions of the PLS landscape, and does an Okay job predicting PLS tree density. Note that the predictions here are overall much higher than the observations for FIA tree denisty.

# Map out predictions in space:

Mapping out predicted - observed for each grid cell in the test dataset.

Plot the predicted - observed by veg class

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Plot predicted - obs against density:

FIA models overpredict over the places that have low tree denisty. I think that if we remove the low denisty places, we might be able to better predict.

Predicting Forest vs. Savanna (not density)

The above models are predicting density (continious variable) from continous envrionmental and climate covariates. However, we can also use the savanna/forest density classificaitons as our y variable and predict the probability a grid cell is forest (ecocode = 1, PLS density > 47 trees/ha) and savannna (ecocode = 0, PLSdensity < 47 trees/ha & >0.5 trees/ha) from continous environmental and climate covariates.

This method has the benefit of predicting the probability of forest (or the probability of savanna) given certain environmental conditions. For this analysis, we exclude prairie sites

Table continues below
  model modeltype formula df AIC
PLS.lgr14 PLS.lgr14 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 50.88 3704
PLS.lgr12 PLS.lgr12 smooth s(sandpct) + s(awc) 27.6 3875
PLS.lgr13 PLS.lgr13 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 27.6 3875
PLS.lgr11 PLS.lgr11 smooth s(sandpct) + s(awc) 35.53 4697
PLS.lgr10 PLS.lgr10 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 27.32 4808
PLS.lgr9 PLS.lgr9 smooth s(sandpct) + s(awc) 27.4 4997
PLS.lgr14L PLS.lgr14L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 5248
PLS.lgr12L PLS.lgr12L linear s(sandpct) + s(awc) 4 5257
PLS.lgr13L PLS.lgr13L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 4 5257
PLS.lgr7 PLS.lgr7 smooth s(sandpct) + s(awc) 18.85 5318
PLS.lgr11L PLS.lgr11L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 5933
PLS.lgr10L PLS.lgr10L linear s(sandpct) + s(awc) 4 6021
PLS.lgr8 PLS.lgr8 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 18.01 6525
PLS.lgr9L PLS.lgr9L linear s(sandpct) + s(awc) 4 6884
PLS.lgr8L PLS.lgr8L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 6888
PLS.lgr5 PLS.lgr5 smooth s(sandpct) + s(awc) 9.597 6894
PLS.lgr5L PLS.lgr5L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 2 6971
PLS.lgr7L PLS.lgr7L linear s(sandpct) + s(awc) 3 7109
PLS.lgr1 PLS.lgr1 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 9.147 7158
PLS.lgr2 PLS.lgr2 smooth s(sandpct) + s(awc) 9.813 7453
PLS.lgr3 PLS.lgr3 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 9.884 7645
PLS.lgr6 PLS.lgr6 smooth s(sandpct) + s(awc) 9.833 7823
PLS.lgr4 PLS.lgr4 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 9.658 7826
PLS.lgr3L PLS.lgr3L linear s(sandpct) + s(awc) 2 7998
PLS.lgr2L PLS.lgr2L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 2 8077
PLS.lgr6L PLS.lgr6L linear s(sandpct) + s(awc) 2 8116
PLS.lgr4L PLS.lgr4L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 2 8148
PLS.lgr1L PLS.lgr1L linear s(sandpct) + s(awc) 2 8414
  MSPE
PLS.lgr14 0.09428
PLS.lgr12 0.09668
PLS.lgr13 0.09668
PLS.lgr11 0.1199
PLS.lgr10 0.1203
PLS.lgr9 0.1246
PLS.lgr14L 0.1346
PLS.lgr12L 0.1352
PLS.lgr13L 0.1352
PLS.lgr7 0.1333
PLS.lgr11L 0.1505
PLS.lgr10L 0.1497
PLS.lgr8 0.1709
PLS.lgr9L 0.1731
PLS.lgr8L 0.18
PLS.lgr5 0.1792
PLS.lgr5L 0.1802
PLS.lgr7L 0.1831
PLS.lgr1 0.1762
PLS.lgr2 0.1915
PLS.lgr3 0.1989
PLS.lgr6 0.1993
PLS.lgr4 0.2018
PLS.lgr3L 0.2092
PLS.lgr2L 0.2095
PLS.lgr6L 0.208
PLS.lgr4L 0.2108
PLS.lgr1L 0.2203

Looks like model number 14 with just Mean Annual Precipitation has lowest AIC for the logistic/binomial model. Lets map out the probability of forest based on model with the lowest AIC value:

forest type ~ s(MAP1910)+ s(pasttmean) + s(deltaT) + s(pastdeltaP)

The smooth predictors in this model are the same as the smooth predictors in the denisty model

savanna forest model for FIA:

  model modeltype formula df AIC MSPE
FIA.lgr14 FIA.lgr14 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 34.44 685.7 0.1092
FIA.lgr12 FIA.lgr12 smooth s(sandpct) + s(awc) 22.58 717.7 0.1066
FIA.lgr13 FIA.lgr13 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 22.58 717.7 0.1066
FIA.lgr11 FIA.lgr11 smooth s(sandpct) + s(awc) 22.15 824.1 0.1366
FIA.lgr10 FIA.lgr10 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 20.14 832.3 0.1367
FIA.lgr9 FIA.lgr9 smooth s(sandpct) + s(awc) 22.53 862 0.1458
FIA.lgr7 FIA.lgr7 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 15.17 888.6 0.1556
FIA.lgr14L FIA.lgr14L linear s(sandpct) + s(awc) 5 1016 0.1397
FIA.lgr11L FIA.lgr11L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 5 1032 0.1537
FIA.lgr12L FIA.lgr12L linear s(sandpct) + s(awc) 4 1035 0.1406
FIA.lgr13L FIA.lgr13L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 4 1035 0.1406
FIA.lgr10L FIA.lgr10L linear s(sandpct) + s(awc) 4 1050 0.1514
FIA.lgr8 FIA.lgr8 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 15.17 1085 0.1993
FIA.lgr2 FIA.lgr2 smooth s(sandpct) + s(awc) 9.557 1088 0.216
FIA.lgr9L FIA.lgr9L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 4 1186 0.1772
FIA.lgr4 FIA.lgr4 smooth s(sandpct) + s(awc) 9.565 1192 0.2224
FIA.lgr8L FIA.lgr8L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 1194 0.1881
FIA.lgr5 FIA.lgr5 smooth s(sandpct) + s(awc) 9.793 1196 0.192
FIA.lgr5L FIA.lgr5L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 2 1204 0.188
FIA.lgr1 FIA.lgr1 smooth s(sandpct) + s(awc) 8.648 1209 0.193
FIA.lgr7L FIA.lgr7L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 3 1234 0.1929
FIA.lgr6 FIA.lgr6 smooth s(sandpct) + s(awc) 9.611 1274 0.229
FIA.lgr3 FIA.lgr3 smooth ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 9.157 1297 0.2314
FIA.lgr2L FIA.lgr2L linear s(sandpct) + s(awc) 2 1323 0.2239
FIA.lgr4L FIA.lgr4L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 2 1329 0.2207
FIA.lgr6L FIA.lgr6L linear s(sandpct) + s(awc) 2 1400 0.2203
FIA.lgr1L FIA.lgr1L linear ecocode ~ s(MAP1910) + s(pasttmean) + s(deltaT) + s(pastdeltaP) + 2 1434 0.2376
FIA.lgr3L FIA.lgr3L linear s(sandpct) + s(awc) 2 1456 0.2319 

## png 
##   2

The above map shows ‘savana or forest’ places where the binomial familiy model predicted between 30-70% probability of forest. There are some places where the model predicts deterministic forest and deterministic savanna.

The main effect on predicted forest here is temperature and temperature seasonality.

It seems like we could make an argument for using the pls density GAMs or the PLS binomial family GAM models and that the analysis of both might end up in supplemental materials.

text left over from previous update

Q2. Does climate explain the distribution of tree density in either the past or the modern landscape?

The past and modern climate space of the Midwest is broad spans a large range of precipitation and temperature.

PLS & FIA tree density are not readily explained by the gradient in Precipitation, but tree denisty bimodality appears to occur between 600-900 mm/year in the past.

PLS & FIA density are not obviously explained as a function of mean annual temperature, but the bimodality in tree density appears to occur between 0-10 DegC in the past:

Q3. Do soil factors/elevation explain tree density in either the past or the modern landscape?

Ans. % Sand, Available Water Content (awc), and saturated hydraulic conductivity (ksat) do not explain tree density overall on their own. Additionally, PC1 and PC2 from a Principle Components Analysis do not explain tree density.

Q4. What is the environmental space where the past was bimodal? How bimodal is the distribution? Where is this located within the N. American Midwest?

Ans. We can look at bins of precipitation to get an idea of where this bimodality coefficient is likely to be highest:

We can look at the bins of temperature to get an idea of where the BC is likely to be highest:

We can also look at the density distribution with bins of the first Principal Component of the environmental data:

We can also look at the density distribution with bins of the second Principal Component of the environmental data:

To fully answer this question of how bimodal the density is, we have to define a bimodality index. There are a few options for this:

  • Dip Test (null = unimodal distribution)
  • Bimodality Coefficient calculated from size, skewness, and kurtosis of the distributions, though there may be issues with both BC and the Dip test (Pfister et al. 2013)
  • BC > 0.556 suggests bimodality
  • BC is calculated with be calculated in the R package modes

Bimodality coefficient is calculated from a distribuiton of data as follows:

BC = \(\frac{(skew^2 + 1)}{(kurtosis + 3)*\frac{(n-1)^2}{(n-2)*(n-3)}}\)

Here we calculated the BC for a bin range of 25mm of precipitation, with non-overlapping and overlapping bins:

We can also calculate BC of the Principal components (this is not as easily interpretable as binning by precipitation):

Where is the bimodality located within the N. American Midwest?

Within the places that were bimodal in the past, and that are bimodal now, we want to know which were savannas and which were forests. We use a classificaiton scheme from Rheumtella et al. (), However there are alternate classificaitons that we could use, such as the scheme Chicago Wilderness has used. We could also make our own classification schemes using kmeans clustering.

Ecotype Rheumtella ChicagoWilderness
Forest > 47 trees/ha > 100 trees/ha
Savanna 0.5-47 trees/ha 10-100 trees/ha
Prairie < 0.5 trees/ha < 10 trees/ha

Ecotypes using the Rheumtella et al. classificaiton

Ecotypes using the Chicago Wilderness classification

Where is the precipitation range where PLS and FIA densities are bimodal?

  • Precipitation climate space where BC > 5.5
  • Using the BC values of non-overlapping bins of width = 25mm
  • Bimodality occurs at both Low precipitation (525-600mm/year) and intermediate precipitation (800-1125 mm/year) in the past
  • Bimodality occurs less frequently on the modern landscape, but also at some low precipitations (575-625 mm/year) and some intemediate precipitations(975 - 1150 mm/year ) on the modern landscape.

Where is the temperature range where PLS and FIA densities are bimodal?

  • Places that are bimodal across temperature ranges (i.e. not linearly predicted by temperature) are mostly located in southern WI, and Indiana and Illinois.
  • This might suggest that some of the bimodality that occurs at intermediate preicipitation could be due to differences in temperature
  • The temperature range where tree density is bimodal is 7-10 degC.

What is the ‘climate space’ in terms of Principal components where PLS and FIA density is bimodal?

  • Temperature and Precipitation dominate loadings of the PC1
  • PC1 shows a similar region of bimodal savanna and forests as the temperature map
  • soil factors (sand, ksat) dominate the loadings of PC2
  • PC2 predicts stable savannas in regions with sandy soils, but the PC2 has bimodalities elsewhere in the midwest
  • When we map the places that are bimodal in both PC1 and PC2, we have an idea of the places that are bimodal after accounting for soils, temperature & precipitation

Principle Component 1 (temperature & precipitation)

Principle Component 2 (soils)

Bimodality in both PC1 and PC2 space (accounting for temp., precip, & soils)

Interpretation:

Within precipiation ranges that have bistablity, differences in temperature may account for differences in tree density. Other regions that are bimodal w.r.t. precipitaiton can be explained by soil factors, such as differences in soil sandiness (PC2). However, a small subset of areas are bimodal in both PC2 and PC1 space.

Summary:

  • PLS density is bimodal, while the modern density is not
  • We may be able to explain some of the distribution in tree density using envrionmental variables (above, and see GAM table)
  • However, none of the environmental variables here fully account for the the bimodality in the PLS data
  • We demonstrate that a large region of the midwest historically could have been either savanna or forest alternative states that are not predicted by precipiation.
  • More of the bimodality in PLS tree distribution is accounted for by tempearature and soils, but not all.
  • The determininants on bimodality
  • Athough historically, this region supported both savanna and forest states, the modern landscape has vitually no bimodal regions. This suggests that fire suppression, land use, and environmental changes that have occurred over the last ~ 150 years may have played a role in the decline of alternative savanna & forest states in the region.

Questions/things to work on:

  • Are there other variables that we should include?
  • Potentially including PDSI from LBDA
  • Figure predicting where the Midwest could be bimodal based on future climate
  • In terms of framing the paper, I need to discuss the regional environmental changes (fire suppresion, land use) in a way that makes it clear that these changes could have removed the endogenous vegetation feedbacks that are important for maintaining alt. stable states.

References